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We present two ready-to-use numerical algorithms to evaluate convex-roof extensions of arbitrary 
pure-state entanglement monotones. Their implementation leaves the user merely with the task 
of calculating derivatives of the respective pure-state measure. We provide numerical tests of the 
algorithms and demonstrate their good convergence properties. We further employ them in order 
to investigate the entanglement in particular few-spins systems at finite temperature. Namely, we 
consider ferromagnetic Heisenberg exchange-coupled spin-| rings subject to an inhomogeneous in- 
plane field geometry obeying full rotational symmetry around the axis perpendicular to the ring 
through its center. We demonstrate that highly entangled states can be obtained in these systems 
at sufficiently low temperatures and by tuning the strength of a magnetic field configuration to an 
optimal value which is identified numerically. 
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I. INTRODUCTION 

Entanglement, one of the most intriguing features of 
quantum mechanics [U [2] , is undoubtedly an indispens- 
able ingredient as a resource to any quantum computa- 
tion or quantum communication scheme The abil- 
ity to (sometimes drastically) outperform classical com- 
putations using multipartite quantum correlations has 
been demonstrated in various theoretical proposals which 
by now have become well known standard examples 
[H [SJ El [7]. Due to the rapid progress in the fields of 
quantum computation, communication, and cryptogra- 
phy, both on the theoretical and the experimental side, 
it has become a necessity to quantify and study the pro- 
duction, manipulation and evolution of entangled states 
theoretically. 

However, this has turned out to be a rather difficult 
task, as the dimension of the state space of a quantum 
system grows exponentially with the number of qudits 
and thus permits the existence of highly nontrivial quan- 
tum correlations between parties. While bipartite en- 
tanglement is rather well understood (see, e.g., [5]), the 
study of multipartite states (with three or more qudits) 
is an active field of research. 

Several different approaches towards the study of en- 
tanglement exist. Bell's original idea [3] that certain 
quantum states can exceed classically strict upper bounds 
on expressions of correlators between measurement out- 
comes of different parties sharing the same state has been 
widely extended and improved to detect entanglement in 
a great variety of states. Entanglement between photons 
persisting over large distances has been demonstrated 
with the use of Bell-type inequalities (see, e.g., Ref. [TU] 
and references therein). Another more recent approach is 
the concept of entanglement witnesses (TU [T^] . These are 
observables whose expectation value is non-negative for 
separable states and negative for some entangled states. 
Thirdly, the concept of entanglement measures is fo- 
cussing more on the quantification of entanglement: if 



state A has lower entanglement than state B, then A 
cannot be converted into B by means of local operations 
and classical communication. Remarkably, there exist in- 
teresting relations between entanglement measures and 
Bell inequalities [T3] on the one hand, and entanglement 
witnesses [TU [11] on the other hand. In this work, we 
focus on the direct evaluation of entanglement measures. 

Among the many features one can demand of such a 
measure, monotonicity is arguably the most important 
one: an entanglement measure should be non-increasing 
under local operations and classical communication (re- 
flecting the fact that it is impossible to create entangle- 
ment in a separable state by these means). A measure 
exhibiting this property is called an entanglement mono- 
tone, with prominent examples being, e.g., the entangle- 
ment of formation [16], the tangle [17], the concurrence 
[T5] or the measure by Meyer and Wallach [T9 . While 
one measure captures certain features of some states es- 
pecially well, other measures focus on different aspects of 
different states. 

Often, entanglement monotones are defined only for 
pure states and are given as analytical expressions of the 
state's components in a standard basis. Unfortunately, 
quantifying mixed-state entanglement is more involved. 
This is somewhat intuitive, since the measure needs to 
be capable of distinguishing quantum from classical cor- 
relations. A manifestation of this difficulty is the fact 
that the problem of determining whether a given density 
matrix is separable or not is apparently very hard and 
has no known general solution for an arbitrary number 
of subsystems with arbitrary dimensions. The ability to 
study mixed-state entanglement is, however, highly de- 
sirable since mixed-states appear naturally due to vari- 
ous coupling mechanisms of the system under examina- 
tion to its environment. There exists a standard way to 
construct a mixed-state entanglement monotone from a 
pure-state monotone, the so-called convex-roof construc- 
tion [20 , but the evaluation of functions obtained in this 
way requires the solution of a rather involved constrained 
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optimization problem (see Sec. 

In this paper, we present two algorithms targeted at 
solving this optimization problem numerically for any 
given convex-roof entanglement measure. In principle, 
these algorithms can also be applied to any optimization 
problem subjected to the same kind of constraints. The 
first algorithm is an extension of a procedure originally 
used to calculate the entanglement of formation [2T1. It 
is a conjugate gradient method exploiting the geometric 
structure of the nonlinear search space emerging from the 
optimization constraint. The second algorithm is based 
on a real parametrization of the search space, which al- 
lows one to carry out the optimization problem in the 
more familiar Euclidean space using standard techniques. 

In the second part of the paper, we use these algo- 
rithms in order to study the entanglement properties of 
a certain type of spin rings. These systems form a gen- 
eralization to N qubits of our previous study, where we 
had only considered the case N — 3 [22] . In the presence 
of an isotropic and ferromagnetic Heisenberg interaction 
and local in-plane magnetic fields obeying a radial sym- 
metry, it can be argued (see Sec. IV and Ref. [22\) that 



the ground state becomes a local unitary equivalent of 
an almost perfect N-partite Greenberger-Horne-Zeilinger 
(GHZ) state [23 



GHZ±) = (|tT...T)±Ui..-i))/^. 



(1) 



Such a system could hence be used for the production 
of highly entangled multipartite states merely by cooling 
it down to low temperatures. One finds, however, that 
the energy splitting between the ground and first excited 
state vanishes in the same limit as the N-partite approx- 
imate GHZ states become perfect, namely for the mag- 
netic field strength going to zero. Therefore, in order to 
quantitatively identify the magnetic field strengths yield- 
ing maximal entanglement at finite temperature, one has 
to study the system in terms of a suitable mixed-state 
entanglement measure. 

The outline of the paper is as follows: In Sec. [H] we 
review how the evaluation of a convex-roof entanglement 
measure is related to a constrained optimization prob- 
lem. We then develop and describe the numerical algo- 
rithms capable of tackling this problem in Sec. HI We 



also present some benchmark tests, comparing our meth- 
ods to another known algorithm. In Sec. |IV| we describe 
the spin rings mentioned earlier and study their entangle- 
ment properties in terms of a convex-roof entanglement 
measure evaluated using our algorithms. We conclude 
our work in Sec. IVl 



II. CONVEX-ROOF ENTANGLEMENT 
MEASURES AS CONSTRAINED 
OPTIMIZATION PROBLEMS 

Given a pure-state entanglement monotone m, the 
most reasonable properties one can demand of a general- 
ization of TO to mixed states are that this generalization 



is itself an entanglement monotone, and that it prop- 
erly reduces to m for pure states. A standard procedure 
which achieves this is the so-called convex-roof construc- 
tion j^ni mj . Given a mixed state p acting on a Hilbert 
space Ti. of finite dimension d, it is defined as 

where 

- {{P^, \i^^)}t=l:S> rank p I C n, 

S S 

P^>0,^P^^^,p^^P^\'lp^){A\j (3) 
i=l i=l 

is the set of all pure-state decompositions of p. Note 
that the pure states ji/'i) are understood to be normal- 
ized. The numerical value of M{p) is hence defined as an 
optimization problem over the set S)(p). 

In order to apply numerical algorithms to this prob- 
lem, 1){p) must be accessible in a parametric way. This 
parametrization is well-known and is often referred to as 
the Schrodinger-H JW theorem [55] , which we briefly 
outline here for the sake of completeness. 

Let St{k,r) denote the set of all fc x r matrices U E 
(pfcxr Tffi^]^ ^]-^g property WU = l^xr, i-e., matrices with 
orthonormal column vectors (hence we have k > r). The 
first part of the Schrodinger-HJW theorem states that 
every U € St{k,r) yields a pure-state decomposition 
fej G S(p) of the density matrix p by the follow- 

ing construction. Let Ai, \xi), * = Ij • • ■ i'' — rankp de- 
note the eigenvalues and corresponding normalized eigen- 
vectors of p, i.e.. 



(4) 



Note that we have Aj > since p is a density matrix and 
as such a positive semi-definite operator. Given a matrix 
U € St{k,r), define the auxiliary states 

r 

li'i) ^^U,jy/Xj\xj), i=l,...,k. (5) 



It is then readily checked that 

1^,) = {i/VpI)\A 



(6) 
(7) 



is indeed a valid decomposition of p into a convex sum of 
k projectors. 

The second part of the theorem states that for any 
given pure-state decomposition {pi, |'(/'i)}i=i of p, there 
exists a [/ € St{k, r) realizing the decomposition by the 
above construction. This guarantees that by searching 
over the set St(k,r) and obtaining the decompositions 
according to the Schrodinger-HJW theorem, we do not 
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'miss out' on any part of the subset of J)(p) with a fixed 
number of states k. The parameterization is thus com- 
plete, i.e., searching the infimum over St{k,r) is equiv- 
alent to searching over all decompositions with fixed so- 
called cardinality k. This allows us to reformulate the 
optimization problem Eq. (|2| as 



M(p)=min inf h(U), 

k>r UeSt(k,r) 



(8) 



where h{U) is the sum on the right-hand side of Eq. ^ 
obtained via the matrix U from p, i.e.. 



Mf/) = ^P.(C/)m(|V'.(C/))). 



(9) 



Note that we have dropped the p-dependence in the 
above expressions, since p is fixed within a particular 
calculation and only the dependence oi h on U is of rel- 
evance in the following. 



It is clear that in a numerical calculation only a fi- 
nite number of different values for k can be investigated. 
However, it is also intuitive to expect that for some large 
enough value of k, increasing the latter even further has 
only marginal effects. In fact, we have observed numer- 
ically that already k — rank p + A yields very accurate 
results in all tests we have performed (also in the ones 
presented in Sec. IIICl, and we have used this choice 



throughout all numerical calculations within this work. 
Note that for a fixed value of k, also all other decompo- 
sitions with cardinality smaller than k are considered as 
well, since the probabilities pi in the elements of 
are allowed to go to zero (with the convention that the 
corresponding states jV'i) are then discarded). 

Since the algorithms presented in the next section will 
both be gradient-based, the derivatives of Eq. ^ with 
respect to the real and imaginary parts of U evaluated 
at U will be required at some point. We state them here 
for the convenience of the reader. They are given by 



dh 



2\iKe{Uki)ni{\MU))) + Y. 

d 

dlmUki ' 2AJm(C/,,M|^fc((7))) + 5] 



dh 



(^) dm 



Re 



Im 



Im 



^R.kl 



■'l.kl 



dm 



91m 

dm 
9 Im 



\i».{U)) 



(10) 
(11) 



r 



where 



III. NUMERICAL ALGORITHMS 



|0fl,fci(C/)) = ^Pk{U)\i\xi) - \i^e{Uki)\MU)) 



(12) 



\<I>im{U)) = \WPk{U)\i\xi) - Xi lmiUki)\MU)) 



(13) 



and superscripts such as in -0^') denote the ith component 
of the state \ip) in an arbitrary but fixed basis. 

As a last remark, we would like to point out that the 
constraint set St{k, r) is, in fact, a closed embedded sub- 
manifold of C'^^'', called the complex Stiefel manifold 
^7J. The geometric structure emerging thereof is ex- 
ploited in one of the two algorithms following shortly. 
The dimension of the Stiefel mainfold is dim St{k,r) = 
2kr — |27l . Since we have k> r,we can set k = r + n, 
n — 0,1, . . .. The number of free parameters N in the op- 
timization is thus N — r^ + 2nr. Hence, N grows linearly 
with n, but quadratically with r. Numerical evaluation 
in larger systems will thus be restricted to low-rank den- 
sity matrices. The fiexibility of choosing n is however 
less restricted. As mentioned above, n = 4 already yields 
satisfying results. 



The study of optimization problems on matrix man- 
ifolds is a rather new and still active field of research 
(see [27j ,28^ and references therein). Only recently, two 
ready-to-use algorithms for minimization over the com- 
plex Stiefel manifold have been presented [221 • To our 
knowledge, these are the only general purpose algorithms 
applicable to generic target functions over St{k, r) found 
in the literature. One is a steepest descent-type method, 
the other one is of Newton-type. We will compare the 
performance of the modified steepest descent algorithm, 
as it is referred to in the original work, with the meth- 
ods presented in this section. We have found that our 
algorithms generally show better convergence properties 
in the cases we have examined. 

We will, however, not make use of the modified New- 
ton algorithm for the following reasons. The second 
derivatives (as required by any Newton-type algorithm) 
of the function h{U) [Eq. (|9|] are in general quite in- 
volved and their number grows quadratically with the 
size of U. Hence, they are very expensive to evaluate, 
even if one resorts to numerical finite differences. More- 
over, the good convergence properties of Newton-type 
methods may only be expected in the very proximity 
of a local minimum. One therefore first typically em- 
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ploys gradient-based techniques to approach a minimum 
sufficiently enough. However, what 'sufficiently enough' 
means in a particular case is often not known before- 
hand. We will later make use of a quasi-Newton algo- 
rithm, which approaches local minima satisfyingly and 
shows strong convergence similar to Newton methods au- 
tomatically when being close enough to a minimum. 



A. Generalized Conjugate-Gradient Method 

fn Ref. [21 a conjugate-gradient algorithm on the uni- 
tary group U{k) = St{k, k) was presented. The goal 
there was to calculate the entanglement of formation also 
for systems with dimensions different from 2x2 [30] . 
Here, we extend this result by noting that the method 
is applicable to any optimization problem on St(k,k), 
particularly to the evaluation of entanglement measures 
other than the entanglement of formation, and we cal- 
culate the required general expression of the gradient of 
h{U). _ ^ 

Optimizing over St{k,k) instead of St(k,r) comes at 
the cost of over-parameterizing the search space. When 
using this algorithm to calculate convex-roof entangle- 
ment measures, we simply took into account only the 
first r columns of the matrix obtained at every iteration. 
This is certainly an aspect one could improve upon in 
future research. 

The algorithm presented here is a conjugate gradient- 
type method, meaning that instead of simply going down- 
hill, i.e., in the direction of steepest descent, previous 
search directions are taken into account at the current 
iteration step. Once the search direction Xi at iteration 
step i, a skew-Hermitian k x k matrix, is known, a line 
search along the geodesic Ui exp{tXi) is performed, where 
Ui is the current iteration point. In particular, one iter- 
ation step of the algorithm may be described as follows 
EH: 



1. Perform a line minimization, i.e., set 



ti+i <- arg min h{U^ exp (tXi)) 



and set 



i+l 



Uiexp{ti+iXi). 



2. Compute the new gradient G^+i at fJ^+i and set 

T ^ exp{t,+iXj2)G^ exp{-t,+iXj2). (16) 

T is the gradient Gi parallel-transported to the new 
point C/i+i. 

3. Calculate the modified Polak-Ribiere parameter 



^ ^ {Gi+i — T, Gi+i) 
{Gi, Gi) 

where {X,Y) = Tr XY^ . 



(17) 



4. Set the new search direction to 

Xi+i < Gi+i + jXi. 



(18) 



5. i 



1. 



6. Repeat from step 1 until convergence. 

The starting point Uq can be chosen arbitrarily, and the 
initial search direction is set to Xq = —Gq. In order to 
find a good approximation to the global minimum, one 
should restart the procedure several times using random 
initial conditions. For the line search in step 1, we uti- 
lized the derivative-free algorithm liiimin described in 
Ref. fST. 

In the following, we calculate the general expression for 
the gradient G of the function h, evaluated at the point 
U (we drop iteration indices for simplicity) . The gradient 
G is defined in terms of the directional derivative of h, 
namely as 



d/i([/(^)(X)) 



de 



{G,X) 



(19) 



e=0 



where U'^'^^X) = V cxp{eX) is a geodesic on St{k,k) 
in direction X (skew-Hermitian matrix) and passing 
through V. The inner product is defined as in step 3 
of the algorithm. We will eventually read off the gradi- 
ent G from its definition in Eq. (19 1. 



Treating h(U) as a function of the real and imaginary 
matrix elements of U, Re Uih and Im Uik, respectively, we 
have 



dh{U^^-^X)) 



de 



£=0 



E 



dh 



dReU, 



ik 



dReU 



dh 
dlmllik 



de 



de 
(e) 

ik 



c=0 



e=Oy 



(20) 



The partial derivatives of h with respect to Re Uik and 



(14) Im [/ifc have already been stated in Eqs. (10 11). Insert 



ing the derivatives of U^^^ into Eq. ( |20| and sorting all 
terms with respect to ReX and ImX, we obtain 



(15) dh{U'^'\X)) 



de 



^Y.^AkiReXki + SkilT^Xki), (21) 



kl 



where 



A 



kl 



E 



dh 



Ski ~ 



dReU,, 



dh 



ReV,k 



dh 



^ V dl-aiUii 



ReK; 



9Im Ui] 



dh 



(22) 



5ReC/,; 



(23) 
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Taking into account the symmetry conditions on X by 
using the relations ReX = {X — X^)/2 and ImX = 
-\{X + X'^) /2 we further obtain 



dh{Ue{X)) 



de 



e=Q 



I E {("^ki - Aik) - i{Ski + Sik))Xki. 



kl 



By comparing this to the right-hand side of Eq. ( 19 ), i.e 
{G,X)^TyGX^ ^-Y^GikXki, 



(24) 
i.e., 

(25) 



kl 



we finally obtain the desired expression for the matrix 
elements of the gradient G, 



Gki = {Aki - Aik)/2 + i{Su + S,k)/2. 



(26) 



One readily sees that G is skew-Hermitian, as required. 

By this, we have completed the description of the 
conjugate gradient algorithm capable of evaluating any 
convex-roof entanglement measure presented in the form 
ofEq. ([8|. 



B. Parametrization with Euler-Hurwitz angles 

Here we present an alternative approach to optimiza- 
tion problems over the Stiefel manifold St{k,r). We will 
obtain a parametrization of St(k, r) in terms of a set 
of real numbers which we will call Euler-Hurwitz an- 
gles, therefore unconstraining the optimization problem 
and mapping it to Euclidean space, where optimization 
problems have been investigated for much longer. We 
will therefore be able to employ a standard algorithm to 
tackle the transformed problem Eq. ([s]) [52] . 

The idea of parameterizing St(k,r) is somewhat moti- 
vated by a theorem known in classical mechanics, where 
it is stated that any rotation in three-dimensional Eu- 
clidean space can be written as a sequence of three el- 
ementary rotations described by three angles, the Euler 
angles. In other words, any orthogonal 3x3 matrix is 
parameterized by three real numbers. It was already Eu- 
ler himself who generalized this idea to arbitrary k x k 
orthogonal matrices i33,J, and Hurwitz [S3] extended the 
parametrization to unitary matrices. We remark that 
ideas in a similar fashion to the ones promoted here 
have been used to calculate an entanglement measure 
for Werner states [^ but were not discussed in greater 
detail. 

We now derive the parametrization of St(k,r). Let 
A G St{k,r). The basic idea is to generate zeroes in 
A and bring it to upper triangular form by applying so- 
called (complex) Givens rotations Gsii^, ifi) [SH] to A from 
the left. The k x k matrices Gs{^, f), s = 1, . . . , fc — 1, 



are defined as 



G^^^{d,ip) 



-e'f sini9, 
e"''^cos??, 
Sij, 



if i ~ j ^ 
if i = s,j = s 
if i = s + 1 , j : 
if i = s + 1, j 
otherwise. 



1, 
s, 
sH 



(27) 



Multiplying A from the left with GsiiSjf), i.e., A = 
Gs{'d,(p)A, has the action 



e'-^cos-d As,. 
-e"^sin -dAs 



e-'fsmi9As 
+ e-''^ cos dA 



if i ~ s, 
if i — s + 1, 
otherwise, 

(28) 



where Ai . denotes the ith row of A. 



Let us write the matrix elements As.j and As+i.j, with 
j arbitrary but fixed, in polar form, i.e., Ag^j = xe"^"' and 
As+i,j = ye"^" , with x, y > 0. We stick to the convention 
that the phases (px and (py be in the interval ]— 7r,7r] in 
order to make this representation unique. It is now easy 
to see that by choosing 



•d — arctan 



0x)/2, 
2/ 

X 



we obtain 



(G,(i9,^)A) 



0, 



(29) 
(30) 

(31) 



while all the other entries in the sth and (s -I- l)th row 
have changed according to Eqs. (28 1. In the case cc = 0, 
we set d = 7r/2 and = Q. In the case j/ = 0, we have 
d and we choose to set = as well. The angles d 
and If are thus restricted to the intervals d € [0, |] and 
LP e]-7r,7r[. 

By successively applying Givens rotations with appro- 
priately chosen angles according to Eqs. (29 1 and (|30|, 



we may now generate zeroes in A column by column, 
from left to right, bottom to top. In greater detail, we 
first erase the whole first column, except for the top en- 
try which will generally remain non-zero. Continuing at 
the bottom of the second column, we may generate zeros 
up to (and including) the third entry from the top of the 
column. If we tried to make the second entry zero, we 
would in general generate a non-zero entry in the second 
row of the first column according to the transformation 
Eq. ( 28 1 . It is convenient to label the angles calculated 



during this process by two indices, and to use the abbre- 
viation Gs{i,j) — Gs{'&ijTfij)- Eventually, we obtain a 
matrix R given by 



R=Q-^A^\{ 



k-l 



n G,{r-i,k-j) 



A. (32) 



The inner of the two products generates zeros in column 
r — i from the bottom up to (and including) row number 
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r — i + 1 . The upper block of R consisting of the first r 
rows is of upper triangular form, while the lower block 
is zero. As a product of unitary Givens rotations, Q^^ 
is itself unitary and in particular invertible. Hence, Q 
always exists and is unitary. We may therefore write 



A=QR= QR, 



(33) 



where Q £ St{k, r) consists of the first r columns of Q 
and R is the upper r x r block of R. Since we assumed 
that A G St{k,r), we have 

l^xr = A^A = {QR)^QR R'^Q'^QR = R^ R, (34) 

and hence, R is unitary. It is straightforward to see that 
a unitary upper triangular matrix can only be of the form 



Rij — ^ij 6 ^ 



(35) 



i.e., a diagonal matrix with only phases on the diagonal. 
Again, we may choose Xi S]"'"'; """l- 

We have thus achieved a unique parametrization of an 
arbitrary matrix A g St{k, r) by a tuple of Euler-Hurwitz 
angles (i?, (/3,x) S S, where 



6 



[0, Ir^'^ 



\r{k- 



(36) 



As required, we find that the number of free parameters 
in this representation is equal to the dimension of the 



Stiefel manifold, i.e., dim St{k,r) — 2kr 



It is clear 



that the procedure described above is fully invertible. 
Hence, we have obtained a one-to-one mapping F : & ^ 
St{k, r). In detail, this mapping, for a vector (t?, ip, x) G 
6, is carried out by filling an otherwise empty kxr matrix 
B with the entries Bu = e'-^% i = 1, . . . r. Then, we apply 
inverse Givens rotations (specified by the Euler-Hurwitz 
angles "i? and tp) from the left to B, in inverse order with 
respect to Eq. (32 1. 



In conclusion, we have transformed the optimization 
problem Eq. (p| into the new problem 



M{p) = min inf h{F{s)). 

k>r sG& 



(37) 



Due to the periodic dependence of F{s) on the angles s, 
it is practical to expand the search space from 6 to the 
whole Euclidean space, making Eq. ( [37| a completely 
unconstrained optimization problem (at the cost of over- 
parameterizing the search space |37jV This problem can 
then be solved using standard numerical techniques. In 
all our calculations, we have used a quasi-Newton algo- 
rithm together with the line search linmin men- 
tioned earlier. This method requires first derivatives 
of the target function with respect to the angles. The 
derivatives with respect to F have already been stated 
in Eqs. (10 11), and the derivatives of F with respect 



to the angles are obtained straightforwardly since each 
angle appears only once in the product representation 
presented above. In order to find a good approximation 
to the global minimum, one should restart with random 
initial conditions several times and take the over-all min- 



C. Test Cases 

Here, we briefly present some performance results of 
the two algorithms presented above. We have applied 
them to the evaluation of two different convex-roof en- 
tanglement measures for which the numerical data can be 
verified by analytically known results. Although our al- 
gorithms show comparatively good performance in these 
cases, we would like to stress that the efficiency of a cer- 
tain method depends strongly on the type of problem 
present, and may even be related to the particular in- 
stance of the problem (see the GHZ /W example below) . 
We have for instance also studied certain matrix approx- 
imation problems, in some of which the parameterized 
quasi-Newton method converged very poorly, whereas 
the modified steepest descent and the generalized con- 
jugate gradient method were equally strong and very ef- 
ficient. One thus cannot generically claim one algorithm 
to be better than the other. It is just beneficial to have 
several different techniques at hand, out of which one 
can choose the best-performing one when applied to a 
particular given problem. 



1. Entanglement of formation of random 2x2 states 

The entanglement of formation |16j is a popular en- 
tanglement measure for bipartite mixed states. It is 
defined as the convex roof of the entropy of entangle- 
ment [3S], which is, for a state the von-Neumann 
entropy S{p) = — Tr plog2 p of the reduced density ma- 
trix p — Ttb Ttb denoting the partial trace over 
the second subsystem. 

Figure [T] shows the convergence behavior of the algo- 
rithms applied to ten random full-rank two-qubit density 
matrices. Displayed is the error at each step of the itera- 
tion between the respective iteration value and the true 
result. The latter is known analytically from Ref. PU] . 

Compared to the algorithms described here, the mod- 
ified steepest descent algorithm due to Ref. [29^ (top 
panel) performs rather poorly. We are aware of the fact 
that we are comparing here a steepest descent algorithm 
with two superlinear algorithms. However, apart from 
presenting convergence properties, we would like to point 
out that the modified steepest descent algorithm often 
converges to imprecise solutions, i.e., it gets stuck in 
undesirable local minima. Rather than on the starting 
point, this phenomenon seems to depend more on the 
actual density matrix itself. 

The conjugate gradient algorithm due to Ref. [21] 
(middle panel) also shows some dependence on the form 
of the density matrix, but always reaches satisfactory ac- 
curacy. The results for the parameterized quasi-Newton 
method (bottom panel) do not, at first glance, show the 
typical fast drop to the solution when close to a good lo- 
cal minimum. This is due to the effect that changing the 
starting point seems to have more influence on the num- 
ber of required iterations in the case of the quasi-Newton 
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FIG. 1: (color online). Convergence plots of the algorithms 
used to evaluate the entanglement of formation on ten ran- 
dom full-rank two-qubit states (each plot was done using the 
same ten states) showing the difference between the numerical 
data and the analytical result as a function of the iteration 
number. Top: The modified steepest descent algorithm from 
Ref. [29] ; Middle: The generalized conjugate gradient method 
from Sec. |III A] Bottom: Quasi-Newton on the parameterized 
search space, Sec. |III B| Each curve in the main plot is av- 
eraged over ten randomly chosen initial points. The typical 
behavior of the algorithms for a single fixed density matrix, 
but with varying initial points of the iteration, is displayed in 
the insets. 



method (see insets in Fig. [T|). When considering single 
(non-averaged) runs of the algorithm, the fast conver- 
gence to the minimum becomes visible. In conclusion, the 
conjugate gradient and the parameterized quasi-Newton 
methods perform best in this case, the latter even slightly 
better than the former. 



2. Tangle of GHZ/W mixtures 

The second test case we present here is concerned with 
the evaluation of the tangle of the rank-2 mixed states 



p(7?) = 7,|GHZ3+)(GHZ+| + (1 - ?7)|W)(W|, 
where |GHZ^) has been defined in Eq. ([T|), 

|w) = -^(im) + im) + iuT)) 



(38) 



(39) 



is the three-qubit W state and < < 1. The tangle 
Tp |T7] is an entanglement measure for pure states of three 
qubits and is known to be an entanglement monotone 
[39] . It can hence be generalized to mixed states by the 
convex roof construction ^ . We will denote the mixed- 
state tangle by t, in contrast to the pure-state version 



T„. The definition of r„ reads 



rp(|^)) = 4 |di-2d2 + 4d3|. 



where 

di 
d2 



(40) 



(41) 

(42) 
(43) 



and ipi,ip2, ■ ■ ■ jips denote the components of the state 
\ip) represented in an arbitrary product basis. In this 
form, the derivatives of Tp with respect to the real and 



imaginary parts of the components of 
the gradient Eqs 



/j), as required by 
(]10 11|), can be read off most easily. 



The tangle takes values between and 1 and is maxi- 
mal for GHZ states. The tangle of the states p{ri) has 
been studied in Ref. [ID], where analytical expressions 
as a function of 77 where presented. Particularly, it was 
found that the tangle vanishes for all < 77 < 770, where 
— 3_^_^3/^ ~ 0.6269, and then continuously increases 



'70 



to unity at rj = 1. 

In Figure [2] we plot the error between the numerically 
obtained and analytically calculated values of T{p{r])) as 
a function of the iteration number for four particular val- 
ues of 7] (see caption of the figure). Only the results of 
the generalized conjugate gradient (top panel) and the 
parameterized quasi-Newton (bottom panel) method are 
shown. The modified steepest descent algorithm from 
Ref. [in] did not succeed to converge to a reasonable local 
minimum for the lowest three values of rj considered. In 
these cases, we empirically find the success rate, which we 
define as the relative number of final errors smaller than 
10~^, tobe<0.1%. For the largest value of 77 examined, 
the algorithm showed typical linear convergence behav- 
ior and arrived at a precision around 10" _ iQ-e ^^j^g^. 
1000 iterations with a rather high success rate of about 
60%. Similarly, the generalized conjugate gradient algo- 
rithm failed to obtain reasonable results for the value of 
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ing this decision. Using the parametrization developed 
in Sec. |IIIB[ one can express each local unitary transfor- 
mation Ui in the matrix U = Ui (S) U2 ^ Un hy its 
Euler-Hurwitz angles and optimize over the whole set of 
all angles simultaneously. Furthermore, on can study in 
this way how 'close' two mixed state are with respect to 
local unitary equivalence. Note that such kind of anal- 
yses are not possible with the modified steepest descent 
or the generalized conjugate gradient methods, since, as 
there is no parametrization, one can optimize over only 
one unitary matrix at a time. 
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FIG. 2; (color online). Convergence of the generalized con- 
jugate gradient (top) and the parameterized quasi-Newton 
(bottom) algorithms for the tangle of GHZ / W states Eq. ( 38 1 
The curves for the values rj = ^rjo (solid line), 77 — (1 — 10 
(dashed line), rj = (1 + 10~*)?7o (dotted line), and rj = 
(dashed-dotted line) have each been obtained by averaging 
100 successful runs starting with random initial points. The 
typical behavior of single runs is shown in the insets, where 
the 100 successful tries that yield one curve in the main plots 
are displayed. 



r] slightly below the threshold value po in most attempts, 
and we find a success rate of < 0.2%. The success proba- 
bility for the other three values of 77 are between 12% and 
95%, whereas they are between 25% and 80% for the pa- 
rameterized quasi-Newton algorithm. One can see, with 
the help of looking more detailed into the behavior of 
single runs (see insets), that the averaged convergence 
plots are slightly fiattened out due to some rather rare 
occurrences of slow convergence. Still, one can observe 
that the parameterized quasi-Newton method converges 
faster to good local minima. 



D. Local unitary equivalence 

We would like to remark here that the parameter- 
ized quasi-Newton method is also capable of determining 
whether two arbitrary mixed states are equivalent up to 
local unitary transformations. While this problem has 
an operational solution in some special cases (see, e.g., 
Ref. im and references therein), there is no generally 
applicable operational criterion known capable of mak- 



IV. PHYSICAL APPLICATION 

In this section, we use the algorithms developed and 
described above to evaluate a multipartite mixed-state 
entanglement measure of a concrete physical system. 



A. Exchange-coupled spin rings with 
inhomogeneous magnetic field geometry 

In the following, we consider the Hamiltonian 



N 



N 



H 



cos ai + Sf sin ai), (44) 



i=l 



where Si = {Sf,Sf,Sf), = 12 with being the 
standard Pauli matrices acting on the ith spin, S^r+i = 
Si and the angles Offe = 27r(fc— 1)/7V, k— 1, . . . , A^. Equa- 
tion (44 1 describes a closed ring of > 2 equidistant 



exchange-coupled spin qubits with local in-plane mag- 
netic fields b,j = (fecos Qfi, &sin ttj, 0)^ which are chosen 
such that the system is invariant under rotations by mul- 
tiples of 2n/N about the center of the ring. The exchange 
coupling J is throughout assumed to be ferromagnetic 
(i.e., J > 0). The fields in Eq. (44 1 are chosen to point 
radially outwards, but the following discussion and re- 
sults also hold for any other local in-plane field configu- 
ration possessing the same rotational symmetry, since all 
these systems are local unitary equivalents. The system 
is depicted schematically in Fig. [s] (a) for three spins. 

In fact, we are considering here a generalization of one 
of the = 3 cases studied in Ref. [32] . There, the partic- 
ular field configuration resulted from semiclassical con- 
siderations with the goal of obtaining a state which is 
close to a GHZ state [see Eq. ([T])] as the ground state of 
the system. In that case, entanglement can be created 
by merely cooling the system to low enough tempera- 
tures. In principle, the argumentation for the occurrence 
of a GHZ ground state presented in Ref. [32] can be ex- 
tended to a number of qubits iV > 3. However, it can 
be expected that for N 00, the lowest-lying multi- 
plet becomes a continuous spectrum. Hence, the ques- 
tion arises up to which numbers of spins TV this setup 
still allows generating GHZ-type entanglement. Before 
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further investigating this question, we briefly restate the 
arguments of Ref. [22] for the convenience of the reader. 

We start from the fact that in the ground state of the 
classical analog of the Hamiltonian (44), all spins are 
aligned for 6=0. However, no direction of alignment 
is favored, reflecting the full rotational symmetry of the 
system in spin space. Small local magnetic flelds (6 <C J), 
applied in the way described above, break this symme- 
try and one is left with the two degenerate ground states 
tt . . . t and ii . . . i where the representation ('quanti- 
zation') axis is the usual z-direction. In fact, each spin is 
slightly tilted against its local magnetic fleld, but there is 
no globally favored direction of orientation, such as with, 
e.g., a global spatially uniform magnetic field. Note that 
this effect of tilting vanishes as 5 — + 0. Due to the Zee- 
man term in Eq. ( 44 1 there is an energy barrier between 



any path connecting the two degenerate minima. In the 
quantum case, tunneling through this barrier lifts the 
degeneracy between the ground states and one obtains 
a tunnel doublet. Thus, in the limit 6^0+, the two 
lowest lying states are the generalized GHZ states given 
in Eq. ([T|. 

As an illustration, we plot the energy surface of the 
classical three-spin system corresponding to Eq. (44 1 in 
Fig.[3](b). We have previously argued [see Ref. [22j, espe- 
cially the discussion leading to Eq. (2) therein] that this 
energy can be expressed in terms of two 'mean' spherical 
ang lesifandd [cf. Fig.[|(a)], since all spins will basically 
align in the present limit 6 <C J, up to small fluctuations 
which sum to zero and are chosen to minimize the total 
energy. One can nicely see how the out-of-plane configu- 
rations at I? = and = vr are energetically favored. For 
any value of ip, a path connecting the two minima has 
to overcome an energy barrier which scales as 0(b^). In 
the figure, this barrier is displayed by the superimposed 
white line for the specific value (p = 'k/2. 

Independently of N, we are generally confronted with 
the following problem if we want to achieve the systems 
considered here to be in a highly entangled state at non- 
zero temperature. On the one hand, the energy splitting 
between the ground state and the first excited state van- 
ishes as b goes to zero. On the other hand, a perfect 
GHZ state is obtained exactly in this limit. For increas- 
ing magnetic field, the states continuously deviate from 
the maximally entangled GHZ state, as can be imagined 
with the help of the classical picture, where the spins 
start to tilt. One therefore has to choose the strength of 
6 as a tradeoff between having a highly entangled ground 
state and separating this state in energy from the next 
higher state. 

In order to find this optimal magnetic field strength at 
a given temperature T ^ we evaluate a suited mixed- 
state entanglement measure on the system's canonical 
density matrix p — exp(— /3iJ)/Trexp(— where /3 — 
l/k^T and fee is Boltzmann's constant. When we stud- 
ied the case iV = 3 in Ref. [22] we used the tangle [see 
Eq. ( [io] )] as our pure-state measure of choice, since it is 
an entanglement measure for three qubits. The general- 
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FIG. 3: (color online), (a) Sch ema tic depiction of the system 
described by the Hamiltonian 1 44 1 for TV = 3. Three spins Si 



are situated at the corners of an equilateral triangle and are 
ferromagnetically exchange-coupled with coupling strength J. 
Local radial in-plane magnetic fields bi (shown as green ar- 
rows in the xy-p\a.ne) point radially outwards. As discussed in 
the text, any other in-plane field geometry obeying the same 
radial symmetry (such as, e.g., a 'chiral' field looping around 
the triangle) leads to equivalent results, (b) Classical energy 
surface Ec of the system shown in the top panel. The 'mean' 
angles '& and (introduced in the top panel) are well suited to 
characterize the state of the system since fluctuations around 
these angles are small for 6 -C J and sum to zero. The super- 
imposed white line shows the perturbatively calculated en- 
ergy barrier at = 7r/2 [see Eq. (2) in Ref. [22]], whereas the 
crosses are due to a corresponding numerical minimization of 
the energy. 



ization to mixed states was done via the convex- roof con- 
struction Eq. ([2]). Here, however, we need a pure-state 
entanglement measure which is defined for any N > 2. 



B. Entanglement measure 

In principle, an exponentially increasing number of dis- 
tinct entanglement measures is required to capture all 
possible quantum correlations in a general pure state of 
N qudits. This may be viewed as the reason for the 
rather large number of proposals for multipartite en- 
tanglement measures that have been put forward over 
the last years. Various insights about the structure and 
characterization of multipartite entanglement have been 
gained by studying such measures. For our purpose, we 
want to have a measure that is easy (and fast) to com- 



10 



pute (in particular, that is an analytic function whose 
complexity grows at most polynomially with N), that 
captures the type of entanglement present in our system 
well, and that possibly has a nice (physical) interpre- 
tation. We found that the Meyer- Wallach measure T9], 
defined for an arbitrary number of qubits, fulfills all these 
criteria. According to Ref. [42 , it can be written in the 
compact form 



7(|V^))=2 



k=l 



(45) 



where pk is the density matrix obtained by tracing out 
all but the fcth qubit out of This is simply the 

subsystem linear entropy averaged over all bipartite par- 
titions involving one qubit and the rest |43j . Moreover, it 
was shown that this entanglement measure is experimen- 
tally observable by determining a set of parameters that 
grows linearly with N, in contrast to the exponentially 
increasing complexity of quantum state tomography f42] . 
We note at this point that the Meyer- Wallach entangle- 
ment has been generalized to a broader family of entan- 
glement measures [44j that might give deeper insight into 
the structure of multipartite entanglement. However, we 
stick to the simple form ( 45 1 for our numerical calcula- 



tions, as this measure turns out to describe our type of 
entanglement well. 

The Meyer- Wallach measure is an entanglement mono- 
tone (and can thus be extended to mixed states via the 
convex- roof construction), lies between zero and one, 
vanishes only for full product states (i.e., states of the 
form |-0) — — 1,...,N), and is maximal for 

gener alized GHZ states Eq. (0. The upper bound is 
however also reached by other states, for instance by 
the so-called cluster states [HI US]. A drawback of the 
Meyer- Wallach measure is that it can also be maximized 
by partially separable states. For example, the state 
1^-) = 1$) (g) 1$), where |$) = (|TT> + \ii))/V2 is a bipar- 
tite Bell state, gives 7(|^')) ~ 1 although it is clearly not 
globally entangled [42 . This is however not a problem in 
our study for two reasons. First of all, we can check by 
numerical diagonalization that the ground state of our 
systems indeed converges to a multiparite GHZ state (at 
least for the first few N < 20). Secondly, comparing the 
data for TV = 3 with our earlier study in Ref. [52] where 
we had employed the tangle, we find the same qualita- 
tive behavior of both entanglement measures. Moreover, 
the optimal values of b for which the measures reach their 
maxima at a given temperature coincide almost perfectly. 
It is thus reasonable to assume that the Meyer- Wallach 
entanglement measure is well suited for quantifying en- 
tanglement in our systems. 

The numerical evaluation of the Meyer- Wallach mea- 
sure extended to mixed states via the convex-roof con- 
struction requires the derivatives of 7(1^/")) with respect 
to the real and imaginary components of [see Eqs. 
(10 11 )]. Due to the partial traces, these expressions are 
a bit cumbersome. However, exploiting the rotational 



symmetry of the Hamiltonian studied here, they can be 
considerably simplified (see Appendix). 



Results 



Before we present and discuss our numerical results, we 



would like to mention that studying the system Eq. ( 44 1 
analytically for arbitrary N is rather difficult. An exact 
diagonalization of the Hamiltonian is not known for ar- 
bitrary N, and perturbation theory to constant order in 
b (independent of N) is not suitable to study the ground- 
state properties of the system, since the ground-state 
splitting is lifted only in A^-th order. One can can thus 
generally expect that the ground-state splitting scales 
with the number of spins as b^ . Since we must always 
have 6^1, this goes to zero for large N, as discussed 
in Sec. |IV A| above. Obtaining highly entangled states 
at finite temperature with this approach will thus be in- 
creasingly difficult for an increasing number of spins N. 

Our numerical results are presented in Figs. [4] and [S] 
Figure |4] shows the Meyer- Wallach measure for N — 
2,3,4, and 5 spins at four different temperatures (see 
caption of the figure). Each data point is the result of 



whichever of the two algorithms described in Sec. HI per- 
formed better in a few trials with random initial condi- 
tions. 

For a fixed number of spins, the entanglement as a 
function of the magnetic field strength b assumes a max- 
imum. This maximal entanglement 7max(?^) is increased 
and its position is shifted to smaller magnetic field val- 
ues as the temperature is lowered. This is due to the fact 
that at low temperatures, only a small magnetic field is 
required in order to make the ground-state splitting suf- 
ficiently large compared with temperature. Since these 
small field values only slightly disturb the ideal GHZ con- 
figuration, almost maximal values of the entanglement 
measure (corresponding to almost perfect GHZ-states) 
are observed. With higher temperature, larger field val- 
ues are required to protect the ground state. Consistent 
with the semiclassical picture, this perturbs the desired 
spin configuration and leads to a lower amount of entan- 
glement. For large magnetic fields, all curves coincide 
eventually, as the system is always found in the ground 
state in that case. 

Figure [s] gives more insight into the dependence of 
maximal entanglement 7max(r) on temperature and the 
number of particles. The plot was obtained by maxi- 
mizing the Meyer- Wallach measure over the magnetic 
field strength b while holding the temperature fixed. Dis- 
played is the difference between the resulting data to the 
zero-temperature maximum (being equal to 1) as a func- 
tion of temperature for different numbers of particles (see 
caption of the figure) . Clearly, the maximally achievable 
entanglement 7max(T') decreases for both increasing tem- 
perature and increasing number of particles. The quali- 
tative dependence on the temperature was discussed al- 
ready above. Here we additionally see an almost linear 
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FIG. 4: (color online). Meyer- Wallach entanglement measure for the system described by the Hamiltonian Eq. (441 at different 
temperatures for several numbers of particles. Concretely, the cases N = 2 (top left), A'^ = 3 (top right). A'' = 4 (bottom 
left), and N = 5 (bottom right) are studied at temperatures T = W~'^J/kB (dashed line), T = IQ-'^J/fca (dashed-dotted line), 
T = IQ-^J/fce (dashed-dot-dotted line), and T = 5 x IQ-^ J/fee (dotted line). 



behavior on a log-log scale at low temperatures, suggest- 
ing a power-law decay of the maximal entanglement of 
the form 1 — 7max(T) oc T" with an exponent a depend- 
ing on the number of particles N . 



The decrease of 7max(T') with the number of spins N at 
fixed temperature is due to the fact that the energy spht- 
ting between the ground and first excited state scales as 

. With a larger number of particles, a higher magnetic 
field is required to achieve a sufficiently large splitting. 
This in turn lowers the entanglement in the ground state, 
due to its 6-dependence, resulting in a lowered maxi- 
mum of the Meyer- Wallach measure. As an additional 
obstacle, the ground-state entanglement as a function of 
b decays even more rapidly as the number of particles 
is increased. This can be seen from the inset of figure |5] 
where, at T = 0, the 6- values yielding the Meyer- Wallach 
measure 0.5 (full width at half maximum, since the max- 
imum at r = is always 1) axe shown as a function of 
iV. 



CONCLUSIONS 



We have presented two ready-to-use numerical algo- 
rithms to evaluate any generic convex-roof entanglement 
measure. While one is based on a conjugate gradient al- 
gorithm operating directly on the search space, the other 
one is a quasi-Newton procedure performing the search 
in the transformed unconstrained Euclidean space. All 
required formulas to implement either of the two algo- 
rithms have been stated explicitly, which, in order to 
calculate different convex-roof extended pure-state mea- 
sures, merely leaves the user with the task of calculat- 
ing its derivatives with respect to the real and imagi- 
nary components of the pure-state argument. The rel- 
atively different nature of the two procedures increases 
the chances that at least one of them performs well in 
the concrete application. In a series of numerical tests, 
we have found that the algorithms perform well and espe- 
cially significantly better than previously presented (non 
Newton-type) ready-to-use optimization problems on the 
Stiefel manifold. However, it is found that the conver- 
gence properties, as is often the case in involved opti- 
mization problems, depend on the cost function. This 
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FIG. 5: (color online). Difference between zero-temperature 
and finite-temperature maximally achievable Meyer- Wallach 
entanglement measure as a function of temperature for sys- 
tems with TV = 2,3,4,5,7,10 spins (from bottom to top). 
Inset: Value of the magnetic field strength 6 as a function 
of system size at which the ground state yields the Meyer- 
Wallach value 0.5 (full width at half maximum). 



N edges. We have argued that a particular local mag- 
netic field geometry, namely radially symmetric in-plane 
fields, favor a highly entangled ground state configura- 
tion. We have confirmed this argumentation by evalu- 
ating the mixed-state Meyer- Wallach entanglement mea- 
sure, defined for an arbitrary number of qubits, and found 
indeed high values of entanglement at low temperatures 
and specific magnetic field strengths. This not only quan- 
tifies the entanglement properties present in this system, 
but also serves more generally as a proof-of-principle for 
the usefulness and applicability of our algorithms. 
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suggests to try applying different techniques to a partic- 
ular optimization problem and examine which one per- 
forms best in that case. 

Further, we have applied our algorithms to evaluate 
a multipartite entanglement measure on density matri- 
ces originating from a real physical system. The latter 
consists of N ferromagnetically exchange-coupled spin- ^ 
particles placed on the edges of a regular polygon with 



APPENDIX: DERIVATIVES OF THE 
MEYER- WALLACH ENTANGLEMENT 
MEASURE 

Within our numerical framework, the evaluation of the 



Meyer- Wallach measure 7(|'0)) [see Eq. (45 1] requires its 
partial derivatives with respect to the real and imaginary 
components of \^p) [see Eqs. ( 10 11 1]. They are given by 



AT 1 1 1 



N 



N 



k—l fi—O 1^1 —0 1^2—0 
AT 1 1 1 1 

EEEE--Ei-(^ 

A;— 1 fi—0 fi—O 1^2—0 



Pk ) 



{v-i_,V2,...,Vk-l,tJ-,Vk-i-l,---:VN) t^k-t^ 

Pk 



(A.l) 
(A.2) 



Here, we represented the zth component of by 



the tuple i 



{vi,...,vm), with the indices Vj G 



{0,1} corresponding to some arbitrary product ba- 
sis {\vi)\v2) ■ ■ ■\vN)}]j-=(i of the spin system. Fur- 
thermore, p^''' denotes the matrix element with in- 
dices (j^, /^) of the reduced density matrix pk = 



Tr. 



Mk-l,Vk + l, 



|^)(^|eC 



2x2 



In case of the systems studied in Sec .|IV A[ the compu- 
tation of the Meyer- Wallach measure (45 1 and its deriva- 



tives can be greatly simplified by exploiting the rota- 
tional symmetry of the Hamiltonian H . Since we have 
\H, TZ\ = 0, where TZ is the symmetry operator for the 
rotation by an angle of 2tt/N about the central axis per- 
pendicular to the plane of the spin ring, all pk in Eq. (45 ) 



are unitary equivalents for a simultaneous eigenstate {ip 
of H and TZ. This reduces Eq. (45 1 to the simple form 



7(|^)) = 2[1-Tr(p?)]. 
The corresponding derivatives read 



(A.3) 



-8[p?'°Re^/.(')-hRe(p;'V 



(2"-l-H^)^ 



(A.4) 



= -8[p?^°ImV'« +Im(p°'V 



(2~-l+^)^ 



(A.5) 
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for < i < 2^-1 - 1, and 
^7 



^7 



1^0) 



IV-) 



(A.6) 
(A.7) 



for 2^^^ < i < 2^ — 1. In practice, we first diagonalize 
H numerically [46 , subsequently diagonalize further any 
degenerate spaces with respect to TZ, and then apply the 
simplified formulas above. 
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